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Local Order and the gapped phase of the Hubbard model: 
plaquette dynamical mean field investigation 
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Abstract. - The four-site "DCA" method of including intersite correlations in the dynamical 
mean field theory is used to investigate the metal-insulator transition in the Hubbard model. 
At half filling a gap-opening transition is found to occur as the interaction strength is increased 
beyond a critical value. The gapped behavior found in the 4-site DCA approximation is shown to 
be associated with the onset of strong antiferromagnetic and singlet correlations and the transition 
is found to be potential energy driven. It is thus more accurately described as a Slater phenomenon 
(induced by strong short ranged order) than as a Mott phenomenon. Doping the gapped phase 
leads to a non-Fermi-liquid state with a Fermi surface only in the nodal regions and a pseudogap 
in the antinodal regions at lower dopings x < 0.15 and to a Fermi liquid phase at higher dopings. 



Understanding the "Mott" or correlation-driven metal 
insulator transition is one of the fundamental questions 
in electronic condensed matter physics [1, 2]. Interest 
increased following P. W. Anderson's proposal that the 
copper oxide based high temperature superconductors are 
doped "Mott insulators" [3]. 1 

Clear theoretical pictures exist in the limits of strong 
and weak coupling. In strong coupling, insulating be- 
havior results from the "jamming" effect [1] in which the 
presence of one electron in a unit cell blocks a second elec- 
tron from entering; we term this the Mott mechanism. At 
weak coupling, insulating behavior arises because long- 
ranged [5] or local [6,7] order opens a gap; we term this 
the Slater mechanism. (Anderson [9] has argued that in 
2d the strong coupling regime provides the appropriate 
description of the low energy behavior for all interaction 
strengths, but this view is controversial and does not ad- 
dress the question of interest here, namely the physical 
origin of the novel low energy physics.) Many materials [2] 



including, perhaps, high temperature superconductors [8] 
seem to be in the intermediate coupling regime in which 
theoretical understanding is incomplete. 

The development of dynamical mean field theory, first 
in its single-site form [10] and subsequently in its cluster 
extensions [11-15] offers a mathematically well-defined ap- 
proach to study metal-insulator transitions. The method, 
while approximate, is non-pcrturbative and provides ac- 
cess to the intermediate coupling regime. In this paper we 
exploit new algorithmic developments [16,17] to obtain de- 
tailed solutions to the dynamical mean field equations for 
the one orbital Hubbard model in two spatial dimensions. 
This, the paradigmatic model for the correlation-driven 
metal-insulator transition, is defined by the Hamiltonian 



H — ^ ^ e p c p,a c P,a + U y ^ 



(1) 



p. a 



It is sometimes useful to distinguish "Mott" materials in which 
the important interaction scale is set directly by an intcrorbital 
Coulomb repulsion from "charge transfer" materials in which the 
interaction scale is set indirectly via the energy required to promote 
a particle to another set of orbitals [4]. For present purposes the 
difference is not important; the term Mott insulator will be used for 
both cases. 



with local repulsion U > 0. We use the electron disper- 
sion e v = —2t(cosp x + cospy). The dynamical mean field 
approximation to this model has been previously consid- 
ered [10,13,18 22]; we comment on the differences to our 
findings below and in the conclusions. 

The dynamical mean field method approximates the 
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electron self energy S(p, uj) by 

E(p,«)= J] 0«(p)£«M. (2) 

a=l...N 

The A functions S a (w) are the self energies of an N- 
site quantum impurity model whose form is specified by 
a self-consistency condition. Different implementations 
of dynamical mean field theory correspond to different 
choices of basis functions <fr a and different self-consistency 
conditions [14, 15]. In this paper we will use primarily 
the "DC A" ansatz [II] although we have also used the 
CDMFT method [12, 22] to verify our results and make 
comparison to other work. In the DCA method one tiles 
the Brillouin zone into N regions, and chooses 4> a (p) = 1 
if p is contained in region a and 4> a (p) = otherwise. 
The "cluster momentum" sectors a correspond roughly to 
averages of the corresponding lattice quantities over the 
momentum regions in which 4> a {p) 0. 

We present results for N = 1 (single-site DMFT) and 
N = 4. Because we are interested in the effects of short 
ranged order, the restriction to small clusters is not a 
crucial limitation: while the precise parameter values at 
which the transition to insulating behavior occurs depend 
on cluster size, the basic relation we establish between cor- 
relations and the insulating behavior does not, and the 4- 
site cluster is computationally manageable so a wide range 
of information can be extracted. 

In the N = 4 case the impurity model is a 4-site cluster 
in which the cluster electron creation operators S may be 
labeled either by a site index j = 1, 2, 3, 4 or by a cluster 
momentum variable A = S, P x , Py, D with S representing 
an average over the range (— 7r/2 < p x < tt/2; —tt/2 < 
p y < 7r/2), P x over the range (n/2 < p x < 3n/2; —n/2 < 
p y < tt/2), and D over the range (tt/2 < p x < 3tt/2; tt/2 < 
Py < 3tt/2). The cluster states are coupled to a bath 
of noninteracting electrons labeled by the same quantum 
numbers. The Hamiltonian is 

H Q I = H cl + J2 (V?<*A, + #bath,(3) 

A, a, a 

A,<r j 

We solve the impurity models on the imaginary frequency 
axis using two new continuous-time methods [16,17]. Be- 
cause we are studying a two dimensional model at temper- 
ature T > we restrict attention to phases without long 
ranged order. The ea, Va an d -f^bath are determined by a 
self consistency condition [10,15]. 

The N = 1 case has been extensively studied [10]. At 
N = 1, intcrsite correlations are entirely neglected; the 
only physics is the strong correlation "local blocking" ef- 
fect envisaged by Mott. If attention is restricted to the 
paramagnetic phase, to temperature T = 0, and density 
n = 1 per site one finds that the ground state is metallic 
for U < U c2 ~ I2t [8] and insulating for U > U c2 . The 
insulating phase is paramagnetic and characterized by an 
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Fig. 1: On-site Green function at time t = (3/2 computed 
using single-site and 4-site DCA methods. All computations 
are performed in the paramagnetic phase at half filling. 

entropy of In 2 per site corresponding to the spin degen- 
eracy of the localized electrons. For U c \ w 9t < U < U C 2 
the insulating phase, although not the ground state, is 
metastable and the extensive entropy of the insulating 
state leads to a transition to the insulating state as the 
temperature is raised [10]. 

The antiferromagnetic solution of the single-site DMFT 
equations has also been extensively studied. The model 
considered here has a nested Fermi surface at carrier con- 
centration n = 1, so at 7i = 1 the ground state is an 
insulating antifcrromagnet at all interaction strengths U. 
The Neel temperature peaks at U « 0.8C/ C 2 [8]. This 
correlation strength also marks a change in the charac- 
ter of the transition: for U < 0.8t/ C 2 the expectation value 
of the interaction term Un^ni decreases as the magnetic 
order increases. The transition is thus potential energy 
driven and is identified with Slater physics. However for 
U > 0.8f/ C 2 the expectation value of the interaction term 
increases as the system enters the antiferromagnetic phase; 
the transition in this case is thus kinetic energy driven and 
is identified with Mott physics. 

We now present results for the N = 4 model in com- 
parison to those obtained in the single-site approxima- 
tion. Figure 1 presents the imaginary time Green function 
G(R,t) at the particular values R = and r = 1/2T = 
(3/2, computed at density n = 1 per site for different tem- 
peratures T and interactions U using 1 and 4 site DCA. 
G(0,/3/2) is directly measured in our simulations and is 
related to the on-site electron spectral function A (lo) by 

The last approximate equality applies for sufficiently small 
T and shows that the behavior of G(0, (3/2) provides infor- 
mation on the existence of a gap in the system. For N = 1 
and U < lOt G(0, (3/2) increases as T decreases, indicating 
the development of a coherent Fermi liquid state. In the 
4-site DCA results a transition is evident as U is increased 
through U* » 4.2i: for U < U* A(0) increases slowly as T 
is decreased, as in the single site model, but for U > U* , 
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Fig. 2: Solid line: on-site spectral function computed by maxi- 
mum entropy analytical continuation of QMC data for U = 6t 
and doping x = 0. Dashed line: spectral function in the 
P — (0, 7r), (it, 0)-momentum sector. Dotted and dash-dotted 
lines: P — (0, tt), (n, 0) and local spectral functions obtained 
by performing the DCA momentum averages of the standard 
SDW mean field expressions for the Green function, with gap 
A = 1.3t. 

A(0) decreases, signaling the opening of a gap. The very 
rapid change across U = U* is consistent with a first or- 
der transition, as found in the careful CDMFT analysis of 
Park et al. [22] . The critical U is seen to be essentially in- 
dependent of temperature indicating that the entropies of 
the metallic and non-metallic states are very similar. The 
end-point of the first order transition is at about T = 0.25t 
which is approximately the Neel temperature of the single- 
site method, at U = At [23]. 

Figure 2 shows as the solid line the local electron spec- 
tral function computed by maximum entropy analytical 
continuation of our QMC data for U = 6i and n = 1. 
Analytical continuation is well known to be an ill-posed 
problem, with very small differences in imaginary time 
data leading in some cases to very large differences in the 
inferred real axis quantities. A measure of the uncertain- 
ties in the present calculation comes from the difference 
between the spectra in the positive energy and negative en- 
ergy regions, which should be equal by particle-hole sym- 
metry. We further note that the gap is consistent with 
the behavior shown in Fig. 1. The local spectral func- 
tion exhibits a characteristic two-peak structure found also 
in CDMFT calculations [22]. The dotted line gives the 
spectral function for the P^-sector, corresponding to an 
average of the physical spectral function over the region 
(7r/2 < p x < 37r/2), (— 7r/2 < p y < 7r/2); this is seen to be 
the origin of the gap-edge structure. 

We present in Fig. 3 the temperature dependence of 
the double-occupancy D = (n^n;) computed using the 
1-site and 4-site DCA for a relatively weak and a rela- 
tively strong correlation strength. In the single-site ap- 
proximation antiferromagnetic correlations are absent in 
the paramagnetic phase and become manifest below the 
Neel temperature; the difference between paramagnetic 
and antiferromagnetic phases therefore gives insight into 



the physics associated with the antiferromagnetic corre- 
lations. For the weaker interaction strength U = 5i, the 
development of Fermi liquid coherence as T is decreased in 
the paramagnetic phase means that the wave function ad- 
justs to optimize the kinetic energy, thereby pushing the 
interaction term farther from its cxtrcmum and increas- 
ing D. At this U the magnetic transition is signaled by 
a rapid decrease in D , indicating that the opening of the 
gap enables a reduction of interaction energy, as expected 
if Slater physics dominates. For the larger U = lOi in the 
single site approximation we see that D is temperature- 
independent in the paramagnetic phase because for this U 
and temperature the model is in the Mott insulating state 
(a first order transition to a metallic state would occur at a 
lower T). The antiferromagnetic transition is signaled by 
an increase in D because in the Mott state the transition 
to antiferromagnetism is kinetic energy driven. 

Turning now to the 4-site calculation we see at U = 5i a 
decrease in D sets in below about T* = 0.23i w 0.8T^ site . 
T* is also the temperature below which G(0,/3/2) begins 
to drop sharply. This indicates that the opening of the gap 
is related to a reduction of interaction energy, implying a 
"Slater" rather than a "Mott" origin for the phenomenon. 
For U = lOt we see a gradual increase in D as T is de- 
creased, reflecting the Mott physics effect of kinetic energy 
gain with increasing local antiferromagnetic correlations. 

To further understand the physics of the transition we 
examine which eigenstates \n c {) oiH c \ are represented with 
high probability in the actual state of the system. We 
define P ncl = (n c i\p c i\n c i) with p c \ the cluster reduced 
density matrix obtained by tracing the partition func- 
tion over the bath states. One particularly interesting 
state is the "plaquette singlet" state which we denote as 
| (12) (34) + (41) (23)) with {ah) representing a singlet bond 
between sites a and b. The upper panel of Fig. 4 shows the 
probability that this state is represented in the thermal en- 
semble corresponding to mean density n = 1 for different 
interaction strengths U; the transition at U « 4.2t mani- 
fests itself as a dramatic change (within our accuracy, the 
jump associated with a first order transition). We have 
performed CDMFT calculations to verify that that the 
same state and same physics control the transition stud- 
ied in Refs. [20,22]. 

The plaquette singlet state has strong intersite corre- 
lations of both d-wave and antiferromagnetic nature. It 
is natural to expect these correlations to open a gap in 
the electronic spectrum. To investigate this possibility 
we computed the DCA momentum averages of the lattice 
Green function using density n = 1, and antiferromagnetic 
and singlet pairing gaps of magnitude A = 1.3t to obtain 
mean field estimates of the impurity model spectral func- 
tions. The dotted and dash-dotted lines in Fig. 2 show 
the antiferromagnetic results. (Use of a <i-wave pairing 
gap would yield very similar results, except that instead 
of a clean gap at one finds a "soft" gap with a linearly 
vanishing density of states). The evident similarity to the 
calculations reinforces the argument that it is the local 
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Fig. 3: Temperature dependence of double occupancy {n-\n{) 
computed using the 1-site and 4-site DCA methods as a func- 
tion of temperature for the half filled Hubbard model at U — 5t 
(upper panel) and U = lOt (lower panel). The 1-site calcula- 
tions are done for both paramagnetic and antiferromagnetic 
phases whereas the 4-site calculation is done for the paramag- 
netic phase only. 
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Fig. 4: Probability that the local Hamiltonian is in a "pla- 
quette singlet state" (a state with plaquette momentum 0) at 
n = 1, T = i/30, as a function of U. The sector statistics are 
measured in the hybridization algorithm. Lower panel: evolu- 
tion of the occupation probabilities with doping at U = 5.2t 
and temperature T = t/30. 



correlations which are responsible for the gapped behav- 
ior. 

We finally consider the effect of doping. The model 
we study is particle-hole symmetric. For definitencss 
we present results for electron doping. In a Fermi liq- 
uid, the imaginary part of the real-axis self energy is 
Im£(p, u> — > 0) oc uj 2 . The spectral representation 
Yi(iu! n ) = J — ImX(p, x)/(iu> n — x) then implies that at 
small uj n , ImE(p, iuj n ) oc oj n . We find that in the S = (0, 0) 
and D = (n, tt) momentum sectors, this relation is obeyed 
at all dopings. The behavior in the P = (0, 7r), (tt, 0)- 
sector is different, as is shown in Fig. 5. The dashed line 
shows the self energy for the half-filled model. The ui^ 1 
divergence, arising from the insulating gap, is evident. For 
large enough doping (x > 0.15) the expected Fermi liquid 
behavior is observed (and indeed for x > 0.2 the self energy 
is essentially the same in all sectors); however for smaller 
dopings, up to x ~ 0.15, ImEp does not extrapolate to 
as uj n — > 0, indicating a non-Fcrmi-liquid behavior in this 
momentum sector. 

To explore further the non-Fermi-liquid behavior we 
present in Fig. 6 the density of states in the P = 
(0, 7r), (tt, 0)-sector, obtained by analytical continuation of 
our quantum Monte Carlo data. Comparison to Fig. 2 



shows that as the chemical potential is increased the Fermi 
level moves into the upper of the two bands. In addition, 
for the lower dopings a small 'pseudogap' (suppression of 
density of states) appears near the Fermi level while for 
x = 0.15 the value of the spectral function at the Fermi 
level approaches that of the noninteracting model, indi- 
cating the restoration of Fermi liquid behavior. We have 
verified that these features are robust, and in particu- 
lar that the suppression of the density of states near the 
Fermi level is required to obtain the measured values of 
G(t ~ /3/2). Comparison of data obtained for inverse 
temperature fit = 30 and j3t = 100 (not shown) with the 
data obtained for fit — 60 shown in Fig. 6 is consistent 
with the pseudogap being the asymptotic low-T behavior, 
not an intermediate T artifact. 

Examination of the D = (tt, 7r)-sector density of states 
and self energy shows that for x = 0.04 and x = 0.08 there 
is no Fermi surface crossing in the D = (tt, 7r)-sector, so 
within the 4-site DCA approximation there is no Fermi 
surface at all. At these chemical potentials most doping 
is provided by incoherent, pseudogapped quasiparticles in 
the P = (0, 7r), (tt, 0)-sector. As x is increased beyond ~ 
0.1 a Fermi crossing appears, first in the D sector and then 
for xj > 0.15 in the P sector, signaling the restoration of 
Fermi liquid behavior. The results may be interpreted as 
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Fig. 5: Imaginary part of Matsubara-axis P — (0, n), (tt, 0)- 
sector self energy measured for U = 5.2t at temperature T = 
t/30 and chemical potential \i (doping x per site) indicated. 




Fig. 6: Doping dependence of P = (0,7r), (ir, 0)-sector density 
of states obtained by analytical continuation of quantum Monte 
Carlo data at U = 5.2t and temperature T = t/60. 



"Fermi arcs" or as hole pockets bounded by the edges of 
the D = (tt, 7r)-sector: the momentum resolution of the 
4-site DCA is insufficient to distinguish the two. As the 
doping is further increased the "Fermi arc" regions rapidly 
grow and the pseudogap fills in, leading to a restoration 
of a conventional Fermi surface for x > 0.15. 

The lower panel of Fig. 4 shows that this non-Fermi- 
liquid behavior can be related to the prominence of the 
plaquette singlet and the plaquette triplet states. The con- 
tribution of the plaquette triplet state peaks at x w 0.15, 
while the contribution of the 6-electron singlet state re- 
mains small, indicating a prominent role for antifcrromag- 
netic (rather than <i-wave singlet) correlations at this dop- 
ing. However, the increasing prominence of the 6-electron 
singlet state as doping is increased strongly suggests that 
the larger doping Fermi-liquid-like state will be suscep- 
tible to a pairing instability. Similar results were found 
in CDMFT calculations by Kyung and collaborators [7], 
who attributed them to antiferromagnetic correlations, by 
Zhang and Imada [20] and by Haule and Kotliar [24]. 

In summary, we have shown that the insulating behav- 
ior (at doping x = 0) and non-Fermi liquid behavior (at 
doping < x < 0.15) found at relatively small U in cluster 
dynamical mean field calculations [7,18-20,22,25] may be 
understood as a consequence of a potential-energy-driven 
transition to a state with definite, strong spatial correla- 
tions, mainly of the plaquctte singlet type. Doping this 
state leads to a low energy pseudogap for momenta in 
the P = (0, 7r), (tt, 0) sector. Superconducting correlations 
(marked by the prominence of the 6 electron states) do not 
become important until beyond the critical concentration 
at which Fermi liquid behavior is restored. Our results 
are consistent with the finding of Park et. al. [22] that the 
JJ-driven transition is first order (although unlike those 
authors we have not performed a detailed study of the 
coexistence region). We interpret the transition as being 
driven by Slater (spatial ordering) physics, whereas Park 
et. al. interpret their results as arising from a strong cou- 
pling, Mott phenomenon. Moukouri and Jarrel [25] argue 



that Slater physics is not important because in a 2d model 
with Heisenberg symmetry long range order does not set 
in until T = 0; We believe, however, that the results for 
double occupancy shown in Fig. 3 and the dominance of 
particular states in the sector statistics plot Fig. 4 provide 
strong evidence that the physics is indeed dominated by 
local order, consistent with Slater-type physics. The im- 
portance of spatial correlations for the spectral function 
and non-Fermi-liquid behavior was previously stressed by 
Jarrell and co-workers [19] and Zhang and Imada [20]. We 
also suggest that the short ranged order is responsible for 
the features noted by Chakraborty and co-workers in the 
optical conductivity and spectral function [21]. Calcula- 
tions in progress will extend the results presented here to 
larger clusters. 
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